% Value Function Iteration for the hybrid menu cost/Calvo model:
%
% Inputs:
%    rPi: The profit function
%    menucost:
%    alphaCalvo: 
%    ErV: An initial guess for ErV (the expected value function)
%    RMU:
%    beta: the discount factor
%    a:
%    rp: The grid for the real price
%    rad:
%    Proba: The probability transition matrix for a (the idiosyncratic shock)
%    ProbNgr:
%    G:
%    agridsize:
%    rpgridsize:
%    radgridsize:
%    tolerance: The error tolerance level of the iteration loop
%    tloop: If one the program prints the time the loop took
%    yesprint: If yesprint == 1 then the function prints 'progress reports'
%
% Outputs:
%    rV: The real valued value function
%    F: The policy function (log price)
%
% Jon Steinsson and Emi Nakamura, August 2006
%********************************************************

function [rV,F,F2] = VFIterationGECalvoPlus(rPi,menucost,...
    menucost2,alphaCalvo,ErV,RMU,beta,a,rp,rad,Proba,...
    ProbNgr,G,tolerance,tloop,yesprint)

agridsize = size(a,1);
rpgridsize = size(rp,1);
radgridsize = size(rad,1);
Ssize = agridsize*rpgridsize*radgridsize;

% GInd is a vector that describes how dividing by inflation moves firms 
% in certain states to new states.
GInd = CalculateGInd(a,rp,rad,G);

GradInd = CalculateGradInd(rad,G);
RMUnew = RMU(GradInd);
RMUnew = reshape(RMUnew,[radgridsize radgridsize]);
ProbNgrNew = ProbNgr.*RMUnew;

time1 = cputime;
i = 1;
ErV_new = ErV;
ErV_old = zeros(agridsize*rpgridsize*radgridsize,1);
Vsubopt = 4*tolerance;
while Vsubopt/2 > tolerance
    ErV_old = ErV_new;
    % Value if nominal price is not changed
    ErV_notchange = rPi + ErV_old;
    % Value if nominal price is changed when menucost is high
    ErV_change = rPi - menucost + ErV_old;
    ErV_change = reshape(ErV_change,[agridsize rpgridsize radgridsize]);
    ErV_change = permute(ErV_change,[2 1 3]);
    ErV_change = reshape(ErV_change,[rpgridsize agridsize*radgridsize]);
    ErV_change = max(ErV_change);
    ErV_change = ones(rpgridsize,1)*ErV_change;
    ErV_change = reshape(ErV_change,[rpgridsize agridsize radgridsize]);
    ErV_change = permute(ErV_change,[2 1 3]);
    ErV_change = reshape(ErV_change,[Ssize 1]);
    % Value if nominal price is changed when menucost is low
    ErV_change2 = rPi - menucost2 + ErV_old;
    ErV_change2 = reshape(ErV_change2,[agridsize rpgridsize radgridsize]);
    ErV_change2 = permute(ErV_change2,[2 1 3]);
    ErV_change2 = reshape(ErV_change2,[rpgridsize agridsize*radgridsize]);
    ErV_change2 = max(ErV_change2);
    ErV_change2 = ones(rpgridsize,1)*ErV_change2;
    ErV_change2 = reshape(ErV_change2,[rpgridsize agridsize radgridsize]);
    ErV_change2 = permute(ErV_change2,[2 1 3]);
    ErV_change2 = reshape(ErV_change2,[Ssize 1]);
    
    % Compare the two and choose the higher when menucost is high
    ErV_mc = (ErV_change > ErV_notchange).*ErV_change ...
        + (ErV_change <= ErV_notchange).*ErV_notchange;
    
    % Compare the two and choose the higher when menucost is low
    ErV_mc2 = (ErV_change2 > ErV_notchange).*ErV_change2 ...
        + (ErV_change2 <= ErV_notchange).*ErV_notchange;
    
    % Take a weighted average of ErV_change and ErV_notchange
    ErV_new = alphaCalvo*ErV_mc + (1-alphaCalvo)*ErV_mc2;
    
    % Take expectations and discount
    ErV_new = (reshape(ErV_new,agridsize,rpgridsize*radgridsize)')*Proba';
    ErV_new = reshape(ErV_new',[Ssize 1]); 
    ErV_new = ErV_new(GInd);
    ErV_new = reshape(ErV_new,[agridsize*rpgridsize radgridsize])*(ProbNgrNew)';
    ErV_new = reshape(ErV_new,[Ssize 1]);
    ErV_new = beta*ErV_new;
    
    % Check tolerance
    Vdifmax = (beta/(1-beta))*max(max(ErV_new - ErV_old));
    Vdifmin = (beta/(1-beta))*min(min(ErV_new - ErV_old));
    Vsubopt = Vdifmax-Vdifmin;
    if Vsubopt/2 <= tolerance
        Vadj = (Vdifmax+Vdifmin)/2;
    end
    
    if yesprint == 1
        if mod(i,5) == 0
            fprintf('.')
        end
        if mod(i,400) == 0
            fprintf('\n')
        end
    end
    i = i + 1;    
end
ErV = ErV_new + Vadj;
clear ErV_new ErV_old ErV_change ErV_change2 ErV_notchange

if tloop == 1
    fprintf('\n\n')
    fprintf('Time of Loop: %8.3f\n\n',cputime - time1)
end

% Create rV
rV_notchange = rPi + beta*ErV;
rp_notchange = repmat(rp,[1 agridsize radgridsize]);
rp_notchange = permute(rp_notchange,[2 1 3]);
rp_notchange = reshape(rp_notchange,[Ssize 1]);

rV_change = rPi - menucost + beta*ErV;
rV_change = reshape(rV_change,[agridsize rpgridsize radgridsize]);
rV_change = permute(rV_change,[2 1 3]);
rV_change = reshape(rV_change,[rpgridsize agridsize*radgridsize]);
[rV_change rp_change] = max(rV_change);
rV_change = ones(rpgridsize,1)*rV_change;
rV_change = reshape(rV_change,[rpgridsize agridsize radgridsize]);
rV_change = permute(rV_change,[2 1 3]);
rV_change = reshape(rV_change,[Ssize 1]);
rp_change = ones(rpgridsize,1)*rp_change;
rp_change = reshape(rp_change,[rpgridsize agridsize radgridsize]);
rp_change = permute(rp_change,[2 1 3]);
rp_change = reshape(rp_change,[Ssize 1]);
rp_change = rp(rp_change,1);

rV_change2 = rPi - menucost2 + beta*ErV;
rV_change2 = reshape(rV_change2,[agridsize rpgridsize radgridsize]);
rV_change2 = permute(rV_change2,[2 1 3]);
rV_change2 = reshape(rV_change2,[rpgridsize agridsize*radgridsize]);
[rV_change2 rp_change2] = max(rV_change2);
rV_change2 = ones(rpgridsize,1)*rV_change2;
rV_change2 = reshape(rV_change2,[rpgridsize agridsize radgridsize]);
rV_change2 = permute(rV_change2,[2 1 3]);
rV_change2 = reshape(rV_change2,[Ssize 1]);
rp_change2 = ones(rpgridsize,1)*rp_change2;
rp_change2 = reshape(rp_change2,[rpgridsize agridsize radgridsize]);
rp_change2 = permute(rp_change2,[2 1 3]);
rp_change2 = reshape(rp_change2,[Ssize 1]);
for ii = 1:length(rp_change2)
    if rp_change2(ii)<=.0
       rp_change2

    end

end

rp_change2 = rp(rp_change2,1);

rV_dif = rV_change - rV_notchange;
rV_dif2 = rV_change2 - rV_notchange;

F_change = (rV_dif>0);
F_change2 = (rV_dif2>0);
F_notchange = (F_change==0);
F_notchange2 = (F_change2==0);
F = F_notchange.*rp_notchange + F_change.*rp_change;
F2 = F_notchange2.*rp_notchange + F_change2.*rp_change2;

rV = alphaCalvo*(F_notchange.*rV_notchange + F_change.*rV_change) ...
    + (1-alphaCalvo)*(F_notchange2.*rV_notchange + F_change2.*rV_change2);
